function S = slcovpca(vmean, C, varargin)
%SLCOVPCA Trains a PCA with the covariance matrix given
%
% [ Syntax ]
%   - S = slcovpca(vmean, C, ...)
%
% [ Arguments ]
%   - vmean:    the mean vector. 
%               It should be either a d x 1 column vector or a zero scalar
%               that indicates a zero mean.
%   - C:        the covariance matrix
%   - S:        the struct of PCA model. (Please refer to slpca for more
%               information about the fields of S)
%
% [ Description ]
%   - S = slcovpca(vmean, C, ...) learns a PCA model from a Gaussian
%     distribution with mean vector and covariance matrix given.
%
%     The following options can be specified to control the learning.
%     \{:
%         - preserve:   The cell array in form of {cri} or {cri, v} to specify
%                       which eigenvalues (and corresponding eigenvectors) are
%                       preserved.
%                       Here, cri is the name of the criteria as listed below, 
%                       v is the parameter value.
%                       \{:
%                           - rank:   preserve all eigenvalues not smaller 
%                                     than eps(max(eigvals))
%                           - num:    preserve the v leading eigenvalues
%                           - ratio:  preserve all eigenvalues not smaller
%                                     than (v * max(eigvals))
%                           - energy: preserve energy of ratio not less
%                                     than v.                           
%                       \:}
%                       By default, preserve is {'rank'}.
%
%         - selfun:     The function handle to select preserved eigenvalues.
%                       It should support the syntax as
%                          $ inds = f(eigvals) $
%                       It outputs the indices of the selected eigenvalues
%                       when input all the eigenvalues.
%                       default = [].
%
%         - maxrank:    The maximum rank of the covariance matrix. 
%                       Typically a covariance matrix computed from n
%                       samples would has a maximum rank at min(d, n-1).
%                       However, n is unknown for such information. 
%
%                       It is recommended to set this option to min(d, n-1)
%                       or simply (n-1) if you know related information in
%                       advance. default = [].
%   
%     \:}
%
% [ Remarks ]
%   - If selfun is specified as a function handle, then the function will
%     use selfun to select the preserved eigenvalues, otherwise the
%     function uses preserve option to select eigenvalues.
%
% [ History ]
%   - Created by Dahua Lin, on Aug 17, 2006
%   - Modified by Dahua Lin, on Jul 17, 2007
%       - totally refactor the program
%       - new options for eigenvalue preservation
%


%% parse and verify input arguments

error(nargchk(2, inf, nargin));

% vmean and C

d = size(C, 1);
assert(isfloat(C) && ndims(C) == 2 && size(C, 2) == d, ...
    'sltoolbox:slcovpca:invalidarg', 'C should be a numeric square matrix.');

assert(isequal(vmean, 0) || (isfloat(vmean) && isequal(size(vmean), [d 1])), ...
    'sltoolbox:slcovpca:invalidarg', 'vmean should be a d x 1 numeric vector.');

% options

opts = struct( ...
    'preserve', {{'rank'}}, ...
    'selfun', [], ...
    'maxrank', []);
if ~isempty(varargin)
    opts = setopts(opts, varargin{:});
end

if isempty(opts.selfun)
    assert(iscell(opts.preserve) && (numel(opts.preserve) == 1 || numel(opts.preserve) == 2), ...
        'sltoolbox:slcovpca:invalidopt', ...
        'The option preserve should be in form of a cell array of length 1 or 2.');    
    fh_seleig = geteigpreserve(opts.preserve{:}); 
else
    assert(isa(opts.selfun, 'function_handle'), 'sltoolbox:slcovpca:invalidopt', ...
        'The option selfun should be either empty or a function handle.');
    fh_seleig = opts.selfun;
end

mrk = opts.maxrank;
if ~isempty(mrk)
    assert(isnumeric(mrk) && isscalar(mrk) && mrk > 0 && mrk == fix(mrk), ...
        'sltoolbox:slcovpca:invalidopt', ...
        'The option maxrank should be either empty or a positive integer scalar.');
end
    

%% main

% compute eigenvalues

[evals, P] = slsymeig(C);

% preprocess

if ~isempty(mrk) && mrk < size(P, 2)
    evals = evals(1:mrk);
    P = P(:, 1:mrk);
end
evals(evals < 0) = 0;
total_energy = sum(evals);

% preserve principal components

inds = fh_seleig(evals);
evals = evals(inds);
P = P(:, inds);

S = makepcamodel(vmean, evals, P, total_energy);

